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Introduction 

We  use  molecular  dynamics  (MD)  both  at  the 
atomistic  and  coarse-grained  level  to  predict  the 
mechanical  and  thermal  properties  of  thermosetting 
polymers.  The  coarse-grained  simulations,  where  the 
polymer  network  is  treated  as  a  bead-spring  system, 
can  capture  several  important  “general”  behaviors  of 
thermosets  such  as  the  role  of  chain  length  of  the 
resin  strands,  degree  of  curing,  strain  rate  and 
temperature  on  the  thermo-mechanical  response  of  a 
cured  polymer  system.  Atomistic  simulations,  on 
the  other  hand,  can  provide  detailed  microscopic 
information  on  the  physical  properties  of 
thermosetting  polymers  and  can  lead  to  predictions 
in  “quantitative”  agreement  with  experiments. 
Recently  a  number  of  MD  models  of  thermosets 
were  developed  and  investigated  by  some 
researchers  [1-7].  While  generating  useful  insight 
into  the  dependence  of  the  physical  properties  of 
thermosetting  polymers  on  their  cross-link  networks, 
these  studies  have  either  not  been  able  to  provide 
specific  correlation  with  the  chemical  structure  of 
the  resin  system  (in  the  case  of  the  coarse-grained 
simulations)  or  been  significantly  affected  by  the 
small  system  size  used  and  the  short  time 
simulations  (in  the  case  of  the  atomistic  studies). 

In  this  proceeding  paper,  we  will  first  briefly  discuss 
our  recent  results,  using  coarse-grained  bead-spring 
model,  on  the  dependence  of  failure  stress  and 
failure  strain  of  highly  cross-linked  polymers  on 
chain  length  of  the  resin  strands.  But,  the  majority  of 
the  paper  will  focus  on  our  work  on  understanding 
the  mechanical  properties  of  one  particular 
thermosetting  polymer,  DGEBA/DETDA  epoxy 
system  using  atomistic  molecular  dynamics 
simulations.  A  series  of  atomistic  simulations  were 


carried  out  to  examine  and  predict  thermo¬ 
mechanical  properties  at  different  extents  of 
reaction.  In  addition,  the  effect  of  crosslink  density 
and  system  size  on  the  properties  of  the  thermosets 
was  also  studied  using  atomistic  simulations. 

Coarse-grained  simulations  of  Thermosets 

We  have  used  a  coarse-grained  bead-spring  model  to 
study  the  dependence  of  the  mechanical  properties 
of  thermosets  on  chain  length  of  the  resin  strands.  In 
the  coarse-grained  model  used  here  the  polymer 
network  is  treated  as  a  bead-spring  system.  To  create 
highly  cross-linked  networks  similar  to  epoxy 
networks,  a  liquid  mixture  consisting  of  strand 
beads,  that  represent  resins,  and  one -bead  molecules, 
representing  cross-linkers,  was  cross-linked 
dynamically.  In  this  investigation  the  strand  length 
was  varied  between  2  and  6  beads  while  the  cross¬ 
linkers  had  a  functionality  of  6.  The  network  was 
dynamically  formed  during  a  constant  temperature 
simulation  until  at  least  95%  of  all  possible  bonds 
between  crosslinkers  and  resin  strands  were  made. 
Details  about  the  coarse-grained  simulation  protocol 
can  be  found  in  [3]. 

The  result  from  coarse-grained  simulation,  shown  in 
Figure  1,  demonstrates  the  relation  between 
mechanical  properties  of  thermosetting  polymers 
and  resin  chain  length.  This  result  clearly  shows  that 
while  the  failure  behavior  of  a  highly  cross-linked 
network  depends  strongly  on  resin  chain  length,  the 
elastic  response  of  the  system  does  not  show  any 
observable  dependence  on  resin  chain  length.  Since 
the  main  objective  of  our  present  work  is  to 
understand  the  elastic  response  of  the 
DGEBA/DETDA  epoxy  system,  the  atomistic 
simulation  part  described  below  will  only  consider  a 
DGEBA  monomer  as  the  resin. 
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For  C  ross-linkers  functionality  of  6 


Figure  1:  Tensile  stress-strain  curves  for  resin 
strands  with  chain  length  2  beads  (black),  3  beads 
(red),  4  beads  (green)  and  6  beads  (blue). 


Building  of  the  Atomistic  models 

The  rest  of  the  paper  will  focus  on  the  atomistic- 
scale  investigation  of  the  thermosetting  materials 
based  on  the  reaction  of  the  DGEBA  (diglycidyl 
ether  of  bisphenol  A)  monomers  with  curing  agent 
EPI-Cure-W  (diethylenetoluenediamine,  DETDA) 
(Figure  2). 


Figure  2:  Reactants  used  for  curing  epoxy-amine 
system,  (left)  Epoxy  resin  monomer  DGEBA 
(diglycidyl  ether  of  bisphenol  A);  (right)  Aromatic 
amine  hardener  DETDA  (diethylene  toluene 
diamine).  The  reaction  sites  are  highlighted  by 
yellow. 


All  the  atomistic  MD  simulations  were  conducted 
using  Discover  modulus  of  the  Materials  Studio 
(Accelrys  Inc.)  software.  The  COMPASS  force  field 
was  employed  in  these  simulations.  The  Velocity 
Verlet  algorithm  was  used  for  integration  of  classical 
Newton’s  equations.  To  build  the  thermosetting 
material,  first,  stoichiometrically  perfect  mixtures  of 
the  epoxy  monomers  and  curing  agent  molecules 
were  packed  into  a  simulation  cell  with  3D  periodic 
boundary  conditions  using  the  Amorphous  Cell 
modulus  of  the  Materials  Studio  package  and  were 
followed  by  energy  minimization  of  the  obtained 
structures.  To  study  the  effect  of  system  size  on  the 


thermo-mechanical  properties  of  the  thermoset 
material  under  consideration,  we  built  mixtures  of 
several  sizes  varying  between  2000  and  8000  atoms. 

After  energy  minimization  of  these  mixtures,  many 
topologically  independent  structures  were  generated 
from  each  of  the  mixtures  by  saving  the 
configuration  at  every  100  ps  during  isothermal  and 
isochoric  (NVT)  equilibration  of  the  mixtures  at 
elevated  temperature  (673  K).  Crosslinking  was  then 
carried  out  on  each  of  the  newly  generated  structures 
(mixtures)  using  the  approach  described  below.  In 
order  to  enhance  molecular  mobility  and  hence  to 
accelerate  the  crosslinking,  the  curing  was  done  at  a 
higher  temperature,  473  K. 

During  the  cross-linking  cycle  chemical  reactions 
occur  based  on  the  distance  criterion.  Bonds  are 
formed  between  reactive  atoms  within  a  cutoff 
distance  and  then  appropriate  hydrogen  atoms  are 
deleted.  To  identify  cross-linking  sites  during  the 
curing  process,  the  capture  sphere  growth  approach 
of  Eichinger  et  al.  [8],  which  allows  specifying  both 
minimum  and  maximum  cutoff  radii  of  the  chemical 
reaction,  was  used.  To  mimic  real  cross-linking 
reactions,  that  is  reactive  groups  have  to  diffuse 
slowly  to  each  other,  the  model  systems  after  each 
cross-linking  cycle  were  subjected  to  a  cascades  of 
both  energy  minimization  and  constant  volume  and 
temperature  MD  equilibrations  according  to  the 
procedure  used  in  Insight  Polymer  software.  The 
polymerization  cycle  was  repeated  until  either  the 
maximum  number  of  cycles  is  reached  or  all  sites 
have  reacted. 

In  order  to  minimize  the  presence  of  defects  and 
stresses  in  a  given  structure,  energetic  and  structural 
information  were  extracted  at  each  cross-linking 
cycle.  Unusually  high  bond  stretching  energy  in  the 
model  system  is  a  sensitive  measure  of  distortion 
and  usually  happens  when  the  cross-linking  process 
introduces  a  defect.  In  the  present  study,  structures 
that  showed  a  dramatic  increase  in  the  maximum 
bond  stretching  energy  at  high  extensions  of  the 
reaction  were  rejected  in  favor  of  stress-free 
structures. 

In  this  work  we  generated  batches  of  many 
topologically  independent  structures  for  each  of  the 
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four  different  system  sizes  studied  (each  containing 
2000,  4000,  6000,  or  8000  atoms).  During  the  curing 
stage  of  each  of  these  structures,  new  structures  with 
conversion  degrees  from  50%  to  100%  were  selected 
and  their  thermo-mechanical  response  was 
characterized. 

Annealing 

As  described  above,  to  enhance  molecular  motions 
and  reduce  stresses  in  the  models,  curing  reaction 
was  conducted  at  the  elevated  temperature  (473  K) 
which  is  above  the  glass  transition  region  for  this 
material.  An  annealing  simulation  procedure  was 
used  to  cool  down  the  systems  with  various  extents 
of  reaction  to  predefined  temperature  values.  In  this 
study  we  chose  two  temperatures  for  mechanical 
characterization  of  the  systems,  298  K  and  480  K 
which  are  above  and  below  the  glass  transition 
temperature  of  a  given  structure,  respectively. 

(i)  Prediction  of  Equilibrium  Properties.  Glass 
Transition  Temperature: 

To  determine  the  glass  transition  temperature  of  a 
given  cross-linked  model,  we  calculated  the  volume- 
temperature  curve  using  constant  pressure  and 
temperature  MD  simulations.  Cooling  of  the  model 
systems  was  carried  out  in  steps  of  10  K  with  100  ps 
equilibration  run  at  each  temperature.  Figure  3 
shows  specific  volume-temperature  plot  that  was 
determined  by  averaging  results  from  5  randomly 
picked  epoxy  structures  that  were  cured  to  the  100% 
extent  of  reaction  and  each  containing  8000  atoms. 


Figure  3:  Specific  volume-temperature  plot  for  fully 
cured  epoxy  structure  containing  8000  atoms. 

The  volume  increases  with  increase  in  temperature 
and  a  change  in  slope  is  observed  for  each  cooling 
cycle.  The  glass  transition  temperatures  were 


determined  at  the  point  of  the  change  in  the  slope 
between  high  and  low  temperature  regions 
(represented  by  the  intersection  of  the  solid  lines  in 
the  above  figure).  The  results  of  our  glass  transition 
investigation  are  presented  in  Figure  4  where  the 
glass  transition  temperature  is  plotted  as  a  function 
of  the  extent  of  the  reaction  for  the  four  different 
system  sizes  studied.  As  expected  [9]  we  observe  an 
increase  in  glass  transition  temperature  with  increase 
in  degree  of  curing.  However,  the  increase  in  Tg  with 
degree  of  curing  is  not  smooth  probably  an 
indication  of  finite  system  size  effect.  The 
experimental  values  of  the  glass  transition 
temperature  for  this  material  found  in  the  literature 
[10-13]  range  from  441  K  to  476  K  and  are  also 
shown  in  figure  4  as  filled  circles.  We  see  that  our 
simulation  results  for  extent  of  reactions  higher  than 
90%  are  within  the  range  of  the  reported 
experimental  values. 
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Figure  4:  Glass  transition  temperatures  as  a  function 
of  the  extent  of  the  reaction.  Triangles  pointing  up: 
2000  atom  structure;  triangles  pointing  down:  4000 
atom  structure;  triangles  pointing  left:  6000  atom 
structure;  triangles  pointing  right:  8000  atom 
structure.  Experimental  values  of  glass  transition 
temperature:  black  filled  circle  [12];  green  circle 

[10]  ;  red  circle  [11] ;  blue  circles  [13]. 

(11)  Coefficent  of  Thermal  Expansion: 

The  slopes  in  the  glassy  and  rubbery  regions  in  the 
specific  volume-temperature  plot  correspond  to  the 
volumetric  coefficients  of  thermal  expansion  (CTE) 
of  the  two  regions.  Figure  5  represents  the  results  on 
coefficients  of  thermal  expansion  calculations  for 
the  glassy  and  rubbery  regions  for  different  degrees 
of  reaction.  In  both  the  glassy  and  rubbery  regions, 
the  CTE  monotonically  decreases  with  the  increase 
in  degree  of  polymerization.  Experimental  value  for 
fully  cured  structure  at  room  temperature  [13]  is 
also  shown  in  figure  5  as  a  black  filled  circle.  The 
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CTE  values  obtained  from  our  simulations  for  high 
degrees  of  curing  are  below  the  experimental  value 
but  are  still  in  a  reasonable  agreement  with 
experiment. 


Figure  5:  Volumetric  coefficients  of  thermal 
expansion  as  a  function  of  the  extent  of  the  reaction 
in  glassy  state  (black  triangles)  and  rubbery  state 
(red  triangles).  The  directions  of  the  triangles 
correspond  to  the  same  model  sizes  as  described  in 
the  Figure  4.  The  filled  circle  represents 
experimental  value  [13]. 

(iii)  Density: 

In  figure  6,  the  dependence  of  density  on  the  extent 
of  curing  reaction  at  temperatures  below  and  above 
Tg  is  shown.  As  expected  the  density  increases  in 
the  course  of  the  reaction  and  the  rate  of  increase 
shows  strong  dependence  on  temperature.  Densities 
of  highly  cured  systems  at  room  temperature  are 
below  but  within  2%  of  the  experimental  value  [14]. 

Mechanical  properties 

Elastic  Coefficents  and  Poisson’s  ratio 

Any  small  volume  element  of  an  amorphous 
material  in  nanoscopic  scale  can  be  characterized  by 
unique  distribution  of  matter  within  it  and  individual 
elements  can  be  viewed  as  regions  of  a 
nanoscopically  heterogeneous  composite  material.  In 
this  sense  we  can  talk  about  distribution  of  the 
material  properties  in  the  macroscopic  sample.  In  the 
present  work  statistical  variation  of  the  elastic 
constants  was  studied  by  estimation  of  upper  and 
lower  bounds  of  the  material  properties  using  Voigt, 
Reuss  and  Hill- Walpole  approaches  [15-18].  For  this 
purpose  we  generated  hutches  of  up  to  60 
topologically  distinct  thermosets  for  each  extent  of 
the  reaction. 
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Figure  6:  Density  as  a  function  of  the  extent  of  the 
reaction  at  298K  (black  triangles)  and  480K  (red 
triangles).  The  directions  of  the  triangles  correspond 
to  the  same  model  sizes  as  described  in  the  Figure  4. 
Filled  circle  represents  experimental  value  [14]. 

To  calculate  elastic  constants  we  used  the  static 
approach  [19].  Three  tensile  and  three  shear 
deformations  of  magnitude  ±0.0005  were  applied  to 
the  systems.  The  obtained  stress  tensor  is  then  used 
to  estimate  stiffness  and  compliance  matrices  Qj  and 
Sy  of  the  material  as  well  as  Young’s,  Shear,  Bulk 
moduli  and  Poisson’s  ratio. 

The  dependence  of  elastic  moduli  on  the  extent  of 
reaction  and  at  two  temperatures  for  the  system 
containing  6000  atoms  is  shown  in  Figure  7.  As 
expected  the  results  from  our  simulation 
demonstrate  the  pronounced  increase  in  Young’s, 
bulk  and  shear  moduli  with  extent  of  curing  at  both 
temperatures.  However,  the  values  of  the  Young’s 
modulus  determined  from  our  simulation  are  above 
the  experimental  values  (2.71  GPa)  even  at  the 
highest  extent  of  the  cross-linking  reaction  [20].  One 
of  the  possible  reasons  of  high  tensile  modulus  of 
the  structures  we  built  for  our  simulations  is  that  our 
models  are  free  of  defects  while  various  kinds  of 
defects  are  usually  present  in  real  macroscopic 
samples  that  tend  to  reduce  the  measured  values  of 
the  elastic  moduli  of  the  material.  As  can  be  seen  in 
figure  7,  the  Poisson’s  ratio  does  not  show  any 
strong  dependence  on  the  extent  of  reaction.  From 
our  simulation  we  determined  an  average  Poisson’s 
ratio  value  of  0.3  for  all  extensions,  which  is  lower 
than  the  experimental  value  [21]. 

There  seems  to  be  a  systematic  dependence  of  the 
elastic  constants  (determined  as  the  average  values 
of  the  upper  and  lower  bounds)  on  system  size, 
where  the  bounds  were  observed  to  narrow 
significantly  with  increase  in  system  size  from  2000 
to  8000  atoms.  We  plan  to  study  large  systems  to 
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minimize  the  effect  of  system  size  so  that  the 
characteristics  of  the  dependence  of  elastic  moduli 
and  other  properties  of  the  thermoset  on  the  extent 
of  the  reaction  can  be  determined  more  precisely. 

Conclusions 

In  the  present  study  we  were  able  to  generate  a  set  of 
stress-free  thermoset  models  with  high  degree  of 
cure.  We  studied  the  dependence  of  the  thermo¬ 
mechanical  properties  of  highly  cross-linked 
polymers  on  the  extent  of  curing  reaction, 
temperature  and  size  of  the  model.  We  found 
densities  and  glass  transition  temperatures  of  the 
systems  in  good  agreement  with  experimental  data. 
However,  the  mechanical  values  from  our  model 
systems  were  found  to  be  higher  than  in  real 
macroscopic  samples.  Size  effect  on  the  material 
properties  was  observed.  Further  increase  in  the  size 
of  the  model  could  help  to  obtain  more  accurate 
simulation  results. 
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Figure  7:  Elastic  moduli  as  a  function  of  the  extent 
of  the  reaction:  (a)  Young's  modulus;  (b)  Poisson ’s 
ratio;  (c)  bulk  modulus;  (d)  shear  modulus.  Black: 
Voigt  and  Reuss  bounds  at  T=298K;  blue:  Hill- 
Walpole  bounds  at  T=298K;  green:  Voigt  and  Reuss 
bounds  at  T=480K;  red:  Hill-Walpole  bounds  at 
T=480K. 
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